# Replication
# John Gerring, Carl Henrik Knutsen, Matthew Maguire, Svend-Erik Skaaning, Jan Teorell, and Michael Coppedge. 2021. “Democracy and Human Development: Issues of Conceptualization and Measurement.” Democratization 28(2): 308-332.

# Clear list objects
rm(list=ls())

# load packages
library(dotwhisker)
library(broom)
library(dplyr)
library(foreign)
library(ggplot2)
library(plm)
library(lmtest)

# Set directory

# Load data
hd_data <- read.dta("2021_Democratization_replication_data_10yr.dta")

# Models

# 1. Polity2 (Marshall)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(e_polity2_pos_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Polity2 (Marshall)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m1.1<-data.frame(term,estimate,std.error,model)
term<-c("Polity2 (Marshall)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m1.2<-data.frame(term,estimate,std.error,model)

# 2. UDS (Pemstein)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(e_uds_mean_pos_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("UDS (Pemstein)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m2.1<-data.frame(term,estimate,std.error,model)
term<-c("UDS (Pemstein)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m2.2<-data.frame(term,estimate,std.error,model)

# 3. Contestation (Miller)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(contdim_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Contestation (Miller)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m3.1<-data.frame(term,estimate,std.error,model)
term<-c("Contestation (Miller)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m3.2<-data.frame(term,estimate,std.error,model)

# 4. Inclusiveness (Miller)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(partdim_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Inclusiveness (Miller)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m4.1<-data.frame(term,estimate,std.error,model)
term<-c("Inclusiveness (Miller)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m4.2<-data.frame(term,estimate,std.error,model)

# 5. Participation (V-Dem)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(v2x_partip_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Participation (V-Dem)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m5.1<-data.frame(term,estimate,std.error,model)
term<-c("Participation (V-Dem)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m5.2<-data.frame(term,estimate,std.error,model)

# 6. Deliberation (V-Dem)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(v2xdl_delib_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Deliberation (V-Dem)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m6.1<-data.frame(term,estimate,std.error,model)
term<-c("Deliberation (V-Dem)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m6.2<-data.frame(term,estimate,std.error,model)

# 7. Female power (V-Dem)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(v2x_gender_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Female power (V-Dem)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m7.1<-data.frame(term,estimate,std.error,model)
term<-c("Female power (V-Dem)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m7.2<-data.frame(term,estimate,std.error,model)

# 8. Civil society (V-Dem)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(v2x_cspart_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Civil society (V-Dem)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m8.1<-data.frame(term,estimate,std.error,model)
term<-c("Civil society (V-Dem)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m8.2<-data.frame(term,estimate,std.error,model)

# 9. Individual liberty (V-Dem)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(v2xcl_rol_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Individual liberty (V-Dem)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m9.1<-data.frame(term,estimate,std.error,model)
term<-c("Individual liberty (V-Dem)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m9.2<-data.frame(term,estimate,std.error,model)

# 10. BMR (Boix)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(e_boix_regime_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("BMR (Boix)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m10.1<-data.frame(term,estimate,std.error,model)
term<-c("BMR (Boix)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m10.2<-data.frame(term,estimate,std.error,model)

# 11. BNR (Bernard)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(e_bnr_dem_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("BNR (Bernard)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m11.1<-data.frame(term,estimate,std.error,model)
term<-c("BNR (Bernard)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m11.2<-data.frame(term,estimate,std.error,model)

# 12. Contestation x Inclusiveness (Miller)
mod.1 <- plm(imr_gapminder_ipo_ln ~ lag(edi_mult_stock_10) + lag(miller_mult_stock_1) + lag(e_migdppcln_ipo), data = hd_data, index = c("country_id","trend"), method = "within", effect = "twoways")
G <- length(unique(hd_data$country_id))
c <- G/(G - 1)
mod.1.cse <- coeftest(mod.1,c * vcovHC(mod.1, type = "HC1", cluster = "group"))
term<-c("Contestation x Inclusiveness (Miller)")
estimate<-c(mod.1.cse[1,1])
std.error<-c(mod.1.cse[1,2])
model <- c("Multiplicative Polyarchy Index (MPI)")
m12.1<-data.frame(term,estimate,std.error,model)
term<-c("Contestation x Inclusiveness (Miller)")
estimate<-c(mod.1.cse[2,1])
std.error<-c(mod.1.cse[2,2])
model <- c("Alternate index")
m12.2<-data.frame(term,estimate,std.error,model)

# plot
dwplot(rbind(m1.1,m1.2,m2.1,m2.2,m3.1,m3.2,m4.1,m4.2,m5.1,m5.2,m6.1,m6.2,m7.1,m7.2,m8.1,m8.2,m9.1,m9.2,m10.1,m10.2,m11.1,m11.2,m12.1,m12.2), dot_args = list(aes(shape = model)), whisker_args = list(aes(linetype = model))) + 
  theme_bw() +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  xlab("Coefficient Estimate") +
  scale_color_grey(start = .1, end = .3) +
  theme(axis.title.x=element_text(size=rel(0.9))) +
  theme(legend.position = "top", legend.title = element_blank(), legend.background = element_rect(size=0.5, linetype="solid",colour ="black")) +
ggsave("Figure_1.pdf")
ggsave("Figure_1.eps")


